Time Series

## Loading required package: s20x

首页

介绍时间序列

时间序列分析：

• 一个时间序列由在连续递增的时间上对数据的收集、重新排序和观察组成。
• 时间间隔相等：时、日、月、季度、周年等
• 时间序列的数据是有序的，不像回归，这里数据的顺序很关键
• 注意，时间是有方向的，如果过去能够影响将来，那么观察值就不是相互独立的。

研究方法

• head(data.df, 5)/tail 查看变量、值的类型
• pairs20x(data.df[,2:5]) 为数据帧绘制直方图（对角线）和散点图（右上角）
• 没有把时间放进去研究的时候
• plot(data.df, which = 1) 查看残差分布，constant and with no trend
• normcheck(lm.fit)
• cooks20x(lm.fit) no influential outlier
• 由于曲线弯曲，所以log是试试
• refit.fit = lm(log(y)~x, data.df)
• summary(lm.fit)
• 查看残差随自变量的分布趋势
• 把时间放进去一起研究的时候残差图看起来更像随机分布在0附近了。可见年顺序在时间序列分析中的作用。

时间序列的数据结构

• 随着时间记录下来的响应变量值
• 以及一系列在同样时间点记录下来的解释变量（表示为：\(X_{t1},X_{t2},X_{t3},...,X_{tk},\))
  • 我们将用它们解释相应变量的行为，或者预测相应变量未来的值。
• 在金融和经济领域，我们通常关心股票市场、消费价格等随时间的变化。
• 这类变量的观察值常常以指数形式出现：
  • 消费价格指数 CPI
  • 交易权重指数 TWI
  • 股指 NZSX50
• 附带解释变量的观察值
  \(时间段 \quad 观察值 \quad 解释变量\)
  \(\quad 1 \qquad Y_1 \qquad X_{11} \quad X_{12} \quad \cdots \quad X_{1k}\)
  \(\quad 2 \qquad Y_2 \qquad X_{21} \quad X_{22} \quad \cdots \quad X_{2k}\)
  \(\quad \vdots \qquad \vdots \qquad \quad \vdots \qquad \vdots \quad \; \vdots \qquad \vdots\)
  \(\quad t \qquad Y_t \; \qquad X_{t1} \; \quad X_{t2} \quad \cdots \quad X_{tk}\)
  \(\quad \vdots \qquad \vdots \qquad \quad \vdots \qquad \vdots \quad \; \vdots \qquad \vdots\)
  \(\quad T \qquad Y_T \qquad X_{T1} \quad X_{T2} \quad \cdots \quad X_{Tk}\)
  
• 这里多了一列时间，它的顺序必须是升序排列

时间序列的成分

时间序列模型

在为时间序列数据建模的时候，分解成四个成分可以帮助理解和量化数据随时间的变化 例如：我们的模型可以如下：
\[Y_t = T_t + C_t + S_T + R_t\]

• 我们可以将这几个成分重新合并成一个新模型来做预测，但这不是理想的方法。

四个成分

1. 趋势成分：
   最慢的成分，我们的问题是：这个趋势是上升的还是下降的？
   • 最慢、长期、平滑
   • 趋势通常来自总体的增长、变化、技术性的改变等等。
   • 趋势可以随着时间或者在一段时间内，上升或下降。
2. 周期成分：
   不规则波浪形行的模式可以出现在不同的时间间隔，并且每次出现都不一样。
   • 第二慢，周期性的，一系列波浪形的
   • 波动的波长和振幅既不固定也无法预测，所以是不规则的模式
   • 很多时间序列的周期性行为主要来自商业周期 注：如果趋势成分允许非线性，则很难区分趋势成分和周期（波动）成分，这样他们就可以合并成一个趋势成分。
3. 季节成分：
   有规律的重复的模式，有固定长度的时间间隔，例如周年或者月。
   • 变化比周期更迅速，为时间序列中，相等间隔规律性的重复模式。
   • 季节性基于日历或者小时
   • 可以字面上理解为三个月（一季度），也可以理解为月，周，日或年。
4. 随机成分(余数)
   无法预测的随机变化
   • 总是存在的，变化最快的成分。
   • 它由那些无法预测的关于趋势，周期，季节成分的变化组成。
   • 但他可以当作其他模型的观测值一样拟合，例如标准正态分布\(N(0, \sigma^2)\), 但它可以不满足独立性原则。
   • 它可以理解为当其他成分都去掉以后留下的成分
   • 它是时间序列中无法被其他东西解释的成分
   • 在季节变化的数据中，随机成分可以从季节成分无法准确重复出现看出。

STL（Seasonal Trend Loess）函数分解时间序列

1. 将响应变量转换成时间序列类型对象 ts()
2. stl()函数执行季节趋势回归分解

data("airpass.df")
passengers.ts = ts(airpass.df$passengers, start = 1949, frequency = 12)
plot(passengers.ts, main = "", ylab = "")

passengers.stl = stl(passengers.ts, s.window = "periodic")
plot(passengers.stl)
abline(v=1950:1960, lty=2)

如何读图

• 解释变量和它的值取自连续的时间间隔
  • 响应变量相对于时间的图有助于发现是否存在时间趋势。

plot(passengers.ts, main = "", ylab = "")

remainder = passengers.stl$time.series[,3]

• • 有对时间进行回归分析得出的残差图可以显示出其他解释变量所没有捕捉到的时间的作用。
  • 将观察值的点连接起来有利于发现时间顺序的模式。plot(y~x, data=data.df) lines(y~x, data=data.df)

例1：

data("mening.df")
mening.ts = ts(mening.df$mening, frequency=12, start=c(1990,1))
plot(mening.ts, main="")
abline(v=1991:2001, lty = 2)

# airpass.df = within(airpass.df, {month = factor(rep(month.abb, 12), levels = month.abb); time = 1:144})
# lines(passengers~month, data = airpass.df)

解读：

• 可以看出有年度季节性模式
• 每隔12个月放一个垂直线，将数据按成年分段：abline(v=1991:2001, lty=2)
• 总体来看这个图是在重复的。
• 所以有季节性。

例2

data(beer.df)
beer.ts = ts(beer.df, frequency = 12, start = c(1971,7))
plot(beer.ts, main = "US Beer Production, July 1970 to June 1978")
abline(v = 1971:1978, lty = 2)

解读：

• 上图总体来看也是重复波动的，而且它的常规峰值出现在一年中相同的时间

例3

passengers.ts = ts(airpass.df$passengers, start = 1949, frequency = 12)
plot(passengers.ts, main = "", ylab="")
abline(v=1950:1960, lty=2)

解读：

• 我们可以看到季节性的涨落幅度更大了；
• 如果这个幅度随时间增加或者减少，我们应该用乘法分解代替加法分解。 \(Y_t = T_t \times C_t \times S_t \times R_t\)
• 由于响应变量是一个计数值，所以用乘法模型是很自然的
• 原则上，我可以通过扩展处理泊松数据的方法，来将其用在计数值时间序列上。
• 这里我们用乘法模型，通过\(log(Y_t)\)来实现，前提是\(Y_t\)都比较大，如本例中。
• 对数据取自然对数可以是季节变化变得恒定，但时间成分的作用将变成倍数（增加多少时间->数据增加多少倍）

log.passengers.ts = log(passengers.ts)
plot(log.passengers.ts, main = "", ylab="")
abline(v=1950:1960, lty=2)

分解图：

• 分解图将帮助我们去掉季节和随机成分，去了解潜在的趋势或波动：

passengers.stl = stl(log.passengers.ts, s.window = "periodic")
plot(passengers.stl)
title("Decomposition of log(Passengers)")

* 趋势存在少量的非线性，最多表明数据中可能有较弱的波动。

分解时间序列：

要解决的问题：

• 相对12月，1月的销量往往有所下降，我们想知道这个下降是否仅仅是正常的季节变化引起的？
• 同样的，由于学生毕业加入劳动力人群，在每个学年末失业率都有所提升。年末提升的失业率表明失业率真的有所提升还是仅仅是正常的季节变化？

解决的方法：

1. 研究去季节化后的值（时间序列去掉季节成分）
   1. 我们估计去季节后的值存在趋势或波动（或者存在解释变量的作用）。研究潜在的趋势或者波动
   2. 如果没有证据证明趋势的存在，我们就可以推断我们看到的只是时间序列上正常的季节变化
   3. 否则我们就说存在根本性的变化，而不是季节
2. 用stl分解
   我们从数据中减掉季节成分再在其上做回归分析；注：这里要小心，因为不满足独立假设。
   1. 获得季节成分
   2. decomp.stl = stl(data.ts, s.window = “periodic”)
   3. sa.data = data.ts - decomp.stl$time.series[,2]

# 代码段2 ----
passengers.ts = ts(log(airpass.df$passengers), start = 1949, frequency = 12)
plot(passengers.ts, main = "Data with Seasonal Component")

passengers.stl = stl(passengers.ts, s.window = "periodic")
sa.passengers = passengers.ts - passengers.stl$time.series[, "seasonal"]
plot(sa.passengers, main = "Data with seasonal component removed")

用一个解释分类变量来去掉季节成分

• 将一个分类变量用在时间序列的回归模型中，分类的级别对应每个季节。
• 用这种方法不需要调整数据去掉季节成分，季节成分是模型设定的一部分
• 这种回归方法适合用在预测，因为它考虑到了整个时间序列。
• 另一种预测方法是Holt-Winters 指数平滑。

预测

指数平滑（exponential smoothing）

概念

1. 指数平滑是能用过去的所有信息来预测将来的方法
2. 每个平滑值是序列中之前所有值的加权平均数（weighted average）
3. 指数序列标记为：\(S_t\) (区别于之前的原始序列\(Y_t\))
4. \(S_t\) 的平滑度有平滑参数\(\alpha\)控制

基本形式：

\[S_t = \alpha Y_t + (1 - \alpha) S_{t-1}\]

• 平滑常量\(\alpha\)越小，结果序列越平滑
• 之所以叫指数平滑是因为距离现在越远的观察值其权重随时间指数级的变小（exponentially smaller）
• 表示如下： \[S_t = \alpha Y_t + \alpha (1 - \alpha) Y_{t-1} + \alpha (1 - \alpha)^2 Y_{t-2} + \cdots + (1-\alpha)^{t-1} Y_1\]
• 从等式的左往右，随着t变得非常大，表达式的系数变得更小（以指数速度 at an exponential rate）

性质

• 平滑后的序列有赖于之前的所有值，其中离现在最近的值的权重最大
• 指数平滑需要大量的观测值
• 以上指数平滑的基本形式不适合存在趋势或者季节性的数据
• 指数平滑在R中用 HoltWinters 函数

指数平滑的基本形式：（废除其他成分，免得成为通用形式）

• HoltWinters(data.df, alpha = x, beta = FALSE, gamma = FALSE) \(\alpha\) 由x来制定，beta 和 gamma 用在有趋势和季节性的时候

data("larain.df")
LArain.ts = ts(rev(larain.df$LA.Rain) * 25.4, frequency = 1, start = 1877)
plot(LArain.ts, xlab = "Season", ylab = "Total rainfall (mm)")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.5, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.5")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.95, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.95")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.75, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.75")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.25, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.25")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.05, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.05")

• 实际上我们可以让HoltWinters来指定\(\alpha\)，它通过最小化均方预测误差来确定。

预测举例：

在predict 中调用 n.ahead 参数

LArain.hw1 = HoltWinters(LArain.ts, beta = FALSE, gamma = FALSE)
LArain.pred = predict(LArain.hw1, n.ahead = 5)
plot(LArain.hw1, predicted.values = LArain.pred, lwd = 2, col.predicted = "blue")

LArain.pred = predict(LArain.hw1, n.ahead = 15, prediction.interval = TRUE)
plot(LArain.hw1, predicted.values = LArain.pred, lwd = 2, col.predicted = "blue")

hist(resid(LArain.hw1))

• 这个模型似乎不是很好，对数据取对数来去掉向右歪斜，同时避免第

log.LArain.ts = log(LArain.ts)
LArain.hw2 = HoltWinters(log.LArain.ts, beta = FALSE, gamma = FALSE)
plot(LArain.hw2)

log.LArain.pred = predict(LArain.hw2, n.ahead = 15, prediction.interval = TRUE)
exp(log.LArain.pred)

## Time Series:
## Start = 1943 
## End = 1957 
## Frequency = 1 
##           fit      upr      lwr
## 1943 343.7816 929.1857 127.1928
## 1944 343.7816 931.5911 126.8644
## 1945 343.7816 933.9966 126.5377
## 1946 343.7816 936.4020 126.2126
## 1947 343.7816 938.8074 125.8893
## 1948 343.7816 941.2129 125.5675
## 1949 343.7816 943.6184 125.2474
## 1950 343.7816 946.0239 124.9289
## 1951 343.7816 948.4296 124.6121
## 1952 343.7816 950.8353 124.2968
## 1953 343.7816 953.2411 123.9831
## 1954 343.7816 955.6470 123.6709
## 1955 343.7816 958.0531 123.3604
## 1956 343.7816 960.4593 123.0513
## 1957 343.7816 962.8657 122.7438

用HoltWinters预测

• 可以用HoltWinters预测包含趋势和季节的时间序列
• 在拟合过程中有三个不同的指数平滑，更新公式：
  • 水平 作用 \(a_t\)
  • 斜率 作用 \(b_t\)
  • 季节 作用 \(s_t\)
• 每一个更新的方程都有它自己平滑常数：\(\alpha, \beta , \gamma\)
• 预测模型为： \[F_{T+h} = a_T + hb_T + s_{T,h}\]
  • 预测是通过估计水平，斜率和季节在时间上的作用来实现的。T是时间序列中的最后一个时间点。
  • \(s_{T,h}\) 表示持续重复的季节作用，尤其是在时间T上估算出来的，并用来计算h个时间单位以后的季节作用。

log.pass.ts = ts(log(airpass.df$passengers), frequency = 12, start = 1949)
log.pass.hw = HoltWinters(log.pass.ts)
log.pass.pred = predict(log.pass.hw, n.ahead = 30)
plot(log.pass.hw, log.pass.pred, col.predicted = "blue", lwd = 2)

让R来决定最优的 \(\alpha, \beta , \gamma\)

log.pass.pred = predict(log.pass.hw, n.ahead = 30, prediction.interval = TRUE)
plot(log.pass.hw, log.pass.pred, col.predicted = "blue", col.intervals = "red", lwd = 2)

针对时间序列的线性回归模型

时间序列线性回归

• 我们前面一直假设在标准模型中，误差是相互独立的，但这个假设在时间序列不成立。
• 通过明确建立依赖模型，可以扩展线性模型来放宽这个假设要求
• 在此简述一阶（first-order）自相关作用

自相关作用

• 拟合线性模型的到时间序列的主要问题就是随机误成分的值常常非独立
• 这个问题成为自相关或者连续相关
• 自相关很难人言识别，尤其是有很多解释变量的时候。所以要集中注意力在解释那些拟合解释变量到响应变量以后所剩下的产值上。

识别自相关

最简单的就是检查残差的自相关函数图（ACF） 首先要明确用回归模型分析时间序列的策略：

1. 以可用的解释变量给均值建模，如果没有，到第二步
2. 通过增加一个季节分类变量来给季节性（如有）建模(这个方法的缺陷是它没有施加任何模式到季节成分上，所以它没有利用1，2月比1月和7月更近这样的事实)
3. 通过增加一个（或多个）时间项到模型中来给所有时间依赖关系建模
4. 检查残差的ACF图。如果有自相关证据，延迟相应（lagged-response）模型可以解决

举例：

plot(log.passengers.ts)

1. 没有解释变量，跳过第一步
2. 有很强的季节效果，建立一个以季节分类建立一个线性模型

# 第二步：添加季节项 ----
airpass.df = within(airpass.df, { month = factor(rep(month.abb, 12), levels = month.abb); time = 1:144})
pass.fit = lm(log(airpass.df$passengers)~month, data = airpass.df)
plot(pass.fit$residuals)
# plot(pass.fit$residuals, type = "l")
lines(airpass.df$time, residuals(pass.fit)) # 用线将点连接起来
lines(lowess(airpass.df$time, residuals(pass.fit))) # 用离散点平滑后画出的曲线，用来了解趋势再好不过了。

# lowess returns a list containing components x and y which give the coordinates of the smooth. The smooth can be added to a plot of the original points with the function lines: see the examples. 在这里显然不用了，因为散点图已经可以看出明显的趋势了。

3. 可以看到很强的时间作用（线性），需要再对时间建模，

# 第三部：添加时间项 ----
pass.fit1 = lm( log(airpass.df$passengers)~time + month, data = airpass.df )
anova(pass.fit1)

## Analysis of Variance Table
## 
## Response: log(airpass.df$passengers)
##            Df  Sum Sq Mean Sq  F value    Pr(>F)    
## time        1 25.1233 25.1233 7143.582 < 2.2e-16 ***
## month      11  2.2843  0.2077   59.047 < 2.2e-16 ***
## Residuals 131  0.4607  0.0035                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-value可以看成指示性的，因为不满足独立性假设。

• 看残差图

plot( residuals(pass.fit1)~fitted(pass.fit1))
lines(lowess(fitted(pass.fit1), residuals(pass.fit1)))

二次方的？

• 重新拟合

pass.fit2 = lm( log(airpass.df$passengers)~time + month + I(time^2), data = airpass.df )
anova(pass.fit2)

## Analysis of Variance Table
## 
## Response: log(airpass.df$passengers)
##            Df  Sum Sq Mean Sq   F value    Pr(>F)    
## time        1 25.1233 25.1233 10813.900 < 2.2e-16 ***
## month      11  2.2843  0.2077    89.386 < 2.2e-16 ***
## I(time^2)   1  0.1587  0.1587    68.307 1.409e-13 ***
## Residuals 130  0.3020  0.0023                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(residuals(pass.fit2)~fitted(pass.fit2))
lines(lowess(fitted(pass.fit2), residuals(pass.fit2)))

plot(residuals(pass.fit2), type = "l") #link the points lower case L
lines(lowess(fitted(pass.fit2), residuals(pass.fit2)))

* 这回看起来好多了，相对恒定而且没有趋势了 * 但残差对时间的图表现出非独立性 4. 残差是否有自相关性？

acf(residuals(pass.fit2))

相关性和自相关性

线性相关系数， \(\rho\)，通常被当做相关性，是衡量两个随机变量\(X 和 Y 之间的线性关系，其中 0<\rho<1\)

• \(如果 \rho = 1 \,则X和Y是完全正相关\)
• \(如果 \rho = -1 则X和Y是完全负相关\)

• \(如果 \rho = 0 \,则X和Y是完全不相关\)

• 我们通过样本相关系数\(r\)来估算相关系数\(\rho\)

• 在简单线性回归中样本的相关系数： \[r = \rm sign(\beta_1) \sqrt{R^2}\]

• 自相关性是变量跟其自身的相关性，但不完全正确，变量Y跟其自身的相关性为1，因为他跟他自己是完全正相关的
• 这里的自相关是延迟h的自相关，是观测值Y和那些h个时间单位以后的Y的值的相关性, T 为时间序列的长度 \[\rm Cor[Y_t, Y_{t+h}]\]
• 例如：h=1时，我们指的相关性是存在于时间t和时间t+1的观测值的相关性

• 这里的相关性是指 \(Y_t 相对于 Y_{t+h}\)的 其中： t=1,…, T-2.

ACF 图

• ACF图是自变量h=0,1,…,T-2，的一个函数图，通常没到T-2就被截断了，因为数据太少了
• 垂线的高度代表我们感兴趣的变量（通常为残差）在时间t和t+h的值的相关性
• 第一条垂直线的高度恒为1，因为\(\rm Cor[Y_t, Y_{t}] = 1\)
• 第二条垂直线的高度是\(\rm Cor[Y_t, Y_{t+1}]\)， lag-1 自相关叫1阶自相关

• 第三条垂直线的高度是\(\rm Cor[Y_t, Y_{t+2}]\)， lag-2自相关又称为二阶自相关，以此类推

• 蓝色的水平虚线是我们认为自相关性I显著的位置。大致来讲，如果垂直线超过虚线，则\(H_0: \rho_h = 0\) 的 P-value将会小于0.05。（由于多重比较，即使没有相关性，有一两个自相关值稍微超过虚线也是正常的）

回到前面如何识别自相关性的第四步

4. 残差是否有自相关性？

acf(residuals(pass.fit2))

* 这是明显的自相关模式，有很多垂直线超过了虚线。我们遇到了自相关问题。 * 要解决这个问题，我们明确建立一个一阶自相关模型，通过引进一个延迟响应变量作为解释变量到模型中来解决这个问题。 * 即：\(Y_1\)用来解释\(Y_2\)，\(Y_2\)用来解释\(Y_3 \cdots\)， 以此类推 * 这意味着我们不能使用\(Y_1\)作为相应变量，因为没有\(Y_0\) log(passengers)[-1] 去掉第一项 * 这类模型称为延迟响应模型

pass.fit3 = lm( log(passengers)[-1] ~ time[-1] + month[-1] + I(time[-1]^2) + log(passengers)[-144], data = airpass.df)

• 我们总要加上延迟相应变量作为模型的最后一项。

anova(pass.fit3)

## Analysis of Variance Table
## 
## Response: log(passengers)[-1]
##                        Df  Sum Sq Mean Sq  F value    Pr(>F)    
## time[-1]                1 24.4515 24.4515 19260.14 < 2.2e-16 ***
## month[-1]              11  2.2733  0.2067   162.79 < 2.2e-16 ***
## I(time[-1]^2)           1  0.1617  0.1617   127.35 < 2.2e-16 ***
## log(passengers)[-144]   1  0.1362  0.1362   107.25 < 2.2e-16 ***
## Residuals             128  0.1625  0.0013                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• P-Value很小，有力证明延迟相应项应该包含在模型中（假设我们的模型是正确的）

acf(residuals(pass.fit3))

• 至少有少量的最高点超过了虚线，这可能是因为多重比较导致的
• 留意一阶自相关建模是如何消除其他显著的自相关的
  • 因为如果A和B紧密相关，而B和C紧密相关，则A和C也可能是紧密相关的。从而消除一阶自相关的时候可能同时消除了更高阶的关系。
  • 有两个延迟（12 和 16）也超过了虚线。这个acf图有效的做了21个假设测试，它给我们图像化展示了前21个自相关性是否显著（在5%水平上）。我们不应该感到奇怪去掉自相关性后还有如果1到两个是显著的。
• 正式的通用延迟响应模型如下： \[Y_t = \beta_0 + \beta_1 t + \beta_2 X_{t1} + \cdots + \beta_{p+1} X_{tp} + \rho Y_{t-1}, \quad t = 2, \cdots , T\]
  • 通过包含一个延迟的响应变量到解释变量中，明确的调整一阶自相关。
  • 延迟相应变量的系数\(\rho\)可以当作残差自相关的天然估值。
  • 这个模型可以用来做预测但不适合作超过往前1步的预测。
• 如果没有修正正自相关就拟合一个线性回归模型，则我们对残差标准差的估值将低于其实际值。
• 结果，将会高估\(R^2\)，还有回归系数的t-test以及anova F-test 也将变得无效。

为什么要关注自相关性

• 信息的遗失：自相关的结果是阻碍我们了解总体中的个体差异的真实大小
• 举例：